/*
 * Software License Agreement (BSD License)
 *
 *  Point Cloud Library (PCL) - www.pointclouds.org
 *  Copyright (c) 2010-2012, Willow Garage, Inc.
 *  Copyright (c) 2012-, Open Perception, Inc.
 *
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions
 *  are met:
 *
 *   * Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *   * Redistributions in binary form must reproduce the above
 *     copyright notice, this list of conditions and the following
 *     disclaimer in the documentation and/or other materials provided
 *     with the distribution.
 *   * Neither the name of the copyright holder(s) nor the names of its
 *     contributors may be used to endorse or promote products derived
 *     from this software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 *  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 *  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 *  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 *  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 *  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 *  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 *  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 *  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 *  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 *
 * $Id$
 *
 */

#ifndef PCL_REGISTRATION_SUPER4PCS_HPP_
#define PCL_REGISTRATION_SUPER4PCS_HPP_

#include <pcl/io/ply_io.h>
#include <pcl/registration/super4pcs.h>
#include <super4pcs/algorithms/super4pcs.h>


struct TransformVisitor {
    inline void operator()(
            float fraction,
            float best_LCP,
            Eigen::Ref<GlobalRegistration::Match4PCSBase::MatrixType> /*transformation*/) const {
      if(fraction >= 0)
        {
          printf("done: %d%c best: %f                  \r",
                 static_cast<int>(fraction * 100), '%', best_LCP);
          fflush(stdout);
        }
    }
    constexpr bool needsGlobalTransformation() const { return false; }
};


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <typename PointSource, typename PointTarget> void
pcl::Super4PCS<PointSource, PointTarget>::computeTransformation (PointCloudSource &output, const Eigen::Matrix4f& guess)
{
  using namespace GlobalRegistration;

  // Initialize results
  final_transformation_ = guess;

  constexpr Utils::LogLevel loglvl = Utils::Verbose;
  using SamplerType   = GlobalRegistration::Sampling::UniformDistSampler;
  using TrVisitorType = typename std::conditional <loglvl==Utils::NoLog,
                            Match4PCSBase::DummyTransformVisitor,
                            TransformVisitor>::type;

  Utils::Logger logger(loglvl);
  MatchSuper4PCS matcher(options_, logger);

  SamplerType sampler;
  TrVisitorType visitor;

  std::vector<GlobalRegistration::Point3D> set1, set2;

  // init Super4PCS point cloud internal structure
  auto fillPointSet = [] (const PointCloudSource& m, std::vector<GlobalRegistration::Point3D>& out) {
      out.clear();
      out.reserve(m.size());

      // TODO: copy other point-wise information, if any
      for(size_t i = 0; i< m.size(); i++){
          const auto& point = m[i];
          out.emplace_back(point.x, point.y, point.z);
      }
  };

  fillPointSet(*target_, set1);
  fillPointSet(*input_, set2);;

  float score = matcher.ComputeTransformation(set1, &set2, final_transformation_, sampler, visitor);
  //if MSSF_Super4PCS algorithm is selected
  if(options_.MSSF){
    std::vector<std::pair<GlobalRegistration::Point3D,GlobalRegistration::Point3D>> congruentSet;
    congruentSet=matcher.getFinalCongruentSet();
    
    for(int i=0;i<4;++i){
      PointSource temPointS;
      PointTarget temPointT;
      temPointT.x=congruentSet[i].first.x();
      temPointT.y=congruentSet[i].first.y();
      temPointT.z=congruentSet[i].first.z();
      base_4Points.push_back(temPointT);
      temPointS.x=congruentSet[i].second.x();
      temPointS.y=congruentSet[i].second.y();
      temPointS.z=congruentSet[i].second.z();
      congruent_4Points.push_back(temPointS);
    }
    //print the congruent set
    cout<<"base_4Points"<<endl;
    ofstream baseKeyPC("base_4Points.txt");
    for(int i=0;i<4;++i){
      baseKeyPC<<base_4Points[i].x<<" "<<base_4Points[i].y<<"  "<<base_4Points[i].z<<endl;
      cout<<base_4Points[i].x<<" "<<base_4Points[i].y<<"  "<<base_4Points[i].z<<endl;
    }
    baseKeyPC.close();
    cout<<"congruent_4Points"<<endl;
    ofstream congruentKeyPC("congruent_4Points.txt");
    for(int i=0;i<4;++i){
      congruentKeyPC<<congruent_4Points[i].x<<" "<<congruent_4Points[i].y<<"  "<<congruent_4Points[i].z<<endl;
      cout<<congruent_4Points[i].x<<" "<<congruent_4Points[i].y<<"  "<<congruent_4Points[i].z<<endl;
    }
    congruentKeyPC.close();

    //MSSF feature calculation
    //obtain the point cloud in R radius
    std::vector<std::vector<int> >pointIdxEightArea;
    pcl::KdTreeFLANN<PointSource>kdtree;
    kdtree.setInputCloud(target_);
    std::vector<int>pointIdxRadiusSearch;
    std::vector<float>pointRadiusSquaredDistance;
    for(int i=0;i<4;i++){
      if(kdtree.radiusSearch(base_4Points[i],options_.spherical_radius,pointIdxRadiusSearch,pointRadiusSquaredDistance)){
        pointIdxEightArea.push_back(pointIdxRadiusSearch);
        pointIdxRadiusSearch.clear();
      }
    }
    kdtree.setInputCloud(input_);
    for(int i=0;i<4;i++){
      if(kdtree.radiusSearch(congruent_4Points[i],options_.spherical_radius,pointIdxRadiusSearch,pointRadiusSquaredDistance)){
        pointIdxEightArea.push_back(pointIdxRadiusSearch);
        pointIdxRadiusSearch.clear();
      }
    }
    //print the PC in radius search
    std::vector<pcl::PointXYZ> radiusSearchPC;
    for(int i=0;i<8;++i){
      if(i<4){
        for(int j=0;j<pointIdxEightArea[i].size();++j){
          pcl::PointXYZ temP=target_->points[pointIdxEightArea[i][j]];
          radiusSearchPC.push_back(temP);
        }
      }
      else{
        for(int j=0;j<pointIdxEightArea[i].size();++j){
          pcl::PointXYZ temP=input_->points[pointIdxEightArea[i][j]];
          radiusSearchPC.push_back(temP);
        }
      }
    }
    ofstream RSPC("radiusSearchPC.txt");
    for(int i=0;i<radiusSearchPC.size();++i){
      RSPC<<radiusSearchPC[i].x<<"  "<<radiusSearchPC[i].y<<" "<<radiusSearchPC[i].z<<endl;
    }
    RSPC.close();

    
    //MSSF feature extraction 
    std::vector< std::vector<std::vector<double> > > MSFList;
    std::vector< std::vector<std::vector<double> > > NorMSFList;
    std::vector< std::vector<Param<Col> > >MSSFList;
    std::vector<pcl::PointIndices> pointIDXList;
    for(int i=0;i<8;i++){
      pcl::PointIndicesPtr ISSIdx(new pcl::PointIndices);
      vector<Param<Col> > atom;
      std::vector<std::vector<double> > FTVector_;
      std::vector<std::vector<double> > NorFTVector_;
      pcl::PointIndices indices;
      indices.indices=pointIdxEightArea[i];
      pcl::PointIndices::Ptr index_Ptr(new pcl::PointIndices(indices));
      if(i<4){
        GenIn::multiScaleSparseFeature(target_, index_Ptr, ISSIdx,0.1,3,FTVector_,NorFTVector_,atom);
      }
      else{
        GenIn::multiScaleSparseFeature(input_, index_Ptr, ISSIdx,0.1,3,FTVector_,NorFTVector_,atom);
      }
      MSFList.push_back(FTVector_);
      NorMSFList.push_back(NorFTVector_);
      MSSFList.push_back(atom);
      pointIDXList.push_back(*ISSIdx);
    }
    //export MSF/NorMSF/MSSF feature
    /**********************************************************/
    for(int i=0;i<MSSFList.size();++i){
      // MSF/NorMSF export
      ofstream outputMSF("MSF"+to_string(i)+".txt");
      ofstream outputNorMSF("NorMSF"+to_string(i)+".txt");
      std::vector<std::vector<double> > temMSF=MSFList[i];
      std::vector<std::vector<double> > temNorMSF=NorMSFList[i];
      for(int j=0;j<temMSF.size();++j){
        for(int k=0;k<temMSF[0].size();++k){
          outputMSF<<temMSF[j][k]<<"  ";
          outputNorMSF<<temNorMSF[j][k]<<"  ";
        }
        outputMSF<<endl;
        outputNorMSF<<endl;
      }
      outputMSF.close();
      outputNorMSF.close();
      //MSSF export
      ofstream output("MSSF"+to_string(i)+".txt");
      vector<Param<Col> > temMSSF;
      temMSSF=MSSFList[i];
      for(int j=0;j<temMSSF.size();++j){
        int rows=temMSSF[0].var.rows();
        for(int k=0;k<rows;++k){
          output<<temMSSF[j].var(k,0)<<"  ";
        }
        output<<endl;
      }
      output.close();
    }
    /**********************************************************/
    //4-points congruent set optimization
    std::vector<std::pair<int,int> > best_matching_pair_list;
    std::vector<std::pair<double,double> > PCCs_pair_list;
    for(int i=0;i<4;++i){
      std::pair<int,int> tem_best_matching_pair;
      std::pair<double,double> tem_PCCs_pair;
      GenIn::mutualMatching(MSSFList[i],MSSFList[i+4],tem_best_matching_pair,tem_PCCs_pair);
      best_matching_pair_list.push_back(tem_best_matching_pair);
      PCCs_pair_list.push_back(tem_PCCs_pair);
    }
    //output the best_matching_pair_list
    ofstream bestMatchingList("best_matching_pair_list.txt");
    for(int i=0;i<best_matching_pair_list.size();++i){
      bestMatchingList<<best_matching_pair_list[i].first<<" "<< 
                      best_matching_pair_list[i].second<<" "<< 
                      PCCs_pair_list[i].first<<" "<< 
                      PCCs_pair_list[i].first<<" "<< endl;
    }
    bestMatchingList.close();
    pcl::registration::TransformationEstimationSVD<pcl::PointNormal,pcl::PointNormal> ::Matrix4 transformation;
    GenIn::transformationCal(target_, input_,best_matching_pair_list,pointIdxEightArea,transformation);
    cout<<"transformation_xu"<<endl<<transformation;
  }
  transformPointCloud (*input_, output, final_transformation_);
  cout<<"final_transformation_: "<<endl<<final_transformation_<<endl;

  pcl::console::print_highlight ("Final score: %f\n", score);

  converged_ = true;
}


#endif

